#Libraries
library(adegenet)
library(vcfR)
library(poppr)
library(ape)
library(viridis)
library(wordcloud)
library(sp)
library(sf)
library(rgdal)
library(adespatial)
library(rgeos)
library(maps)
library(maptools)
library(ggplot2)
library(pals)
library(rcartocolor)
library(dbscan)
library(plotly)
library(gg3D)
library(akima)
library(vegan)
library(hierfstat)
library(pegas)
library(dplyr)
library(tidyr)
library(raster)
library(groc)
library(qvalue)
library(fields)
library(gridExtra)
library(magick)
library(rasterVis)
library(pdftools)
library(ecodist)#Once genotype data is read into R is must be converted into a genind object or other objects (e.g. genlight or genpop) to be useable by adegenet, adespatial, and poppr
#reading in vcf file
tansy_vcf <- read.vcfR("data/stacks_populations_output/populations.snps.vcf")## Scanning file to determine attributes.
## File attributes:
## meta lines: 14
## header_line: 15
## variant count: 3071
## column count: 187
##
Meta line 14 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
## Character matrix gt rows: 3071
## Character matrix gt cols: 187
## skip: 0
## nrows: 3071
## row_num: 0
##
Processed variant 1000
Processed variant 2000
Processed variant 3000
Processed variant: 3071
## All variants processed
pops <- read.csv("data/tansy_pop_info_allsamples.csv")
#converting vcf data object into a genind object that is used in adegent and adespatial
tansy_genind <- vcfR2genind(tansy_vcf)
pop(tansy_genind) <- pops[,2]
ploidy(tansy_genind) <- 2
rownames(tansy_genind@tab) <- pops[,3]
#Fill NAs with mean allele frequency
tansy_tab <- tab(tansy_genind, NA.method = "mean")
#convert into other data objects that are used by poppr and adegenet/adespatial
tansy_genclone <- poppr::as.genclone(tansy_genind)
tansy_genlight <- vcfR2genlight(tansy_vcf)
tansy_loci <- vcfR2loci(tansy_vcf)
tansy_genpop <- genind2genpop(tansy_genind, quiet = TRUE)
tansy_hierfstat <- genind2hierfstat(tansy_genind)#Plot of sites (mostly MN)
#Read in coordinates
xy <- pops[,5:6] #lat/lon of each of the populations in the dataset
coordinates(xy) <- ~ lon + lat
proj4string(xy) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")
# plot(xy, xlim = c(-93, -92), ylim = c(43, 50), pch = 19, cex = 0.1)
# textplot(xy@coords[,1], xy@coords[,2], words = rownames(tansy_tab), cex = 0.25, new = FALSE)
pops$geo_grp <- as.factor(pops$geo_grp)
pops$ECS_fac2 <- as.factor(pops$ECS_fac2)
pops_df <- pops
coordinates(pops) <- ~ lon + lat
proj4string(pops) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")#For use in population level stats where Finland is removed
tansy_genind_pop <- tansy_genind
tansy_genind_pop$other$xy <- xy@coords
nums <- seq(1:178)
fin_pops <- grep("FINLAND", pops_df$name)
wis_pops <- grep("WIS", pops_df$name)
nonMN_pops <- c(fin_pops, wis_pops)
MN_which <- which(!is.element(nums, nonMN_pops))
MN_pops <-MN_which
pops_df <- pops_df[MN_pops,]
xy_MN <- as_tibble(xy@coords[MN_pops,])
#write.csv(xy_MN, file = "data/xy_MN.csv", row.names = FALSE)
xy_MN$ID <- MN_pops
coordinates(xy_MN) <- ~ lon + lat
proj4string(xy_MN) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")
tansy_genind_pop_MN <- tansy_genind_pop[MN_pops,]
tansy_genpop_pop_MN <- genind2genpop(tansy_genind_pop_MN, quiet = TRUE)
tansy_hierfstat_pop_MN <- genind2hierfstat(tansy_genind_pop_MN)midwest <- map_data("county", "minnesota")
midwest <- midwest[midwest$long > -100 & midwest$long < -89.5 &
midwest$lat > 41 & midwest$lat < 52,]
MNsites <- ggplot() + geom_path(data = midwest, mapping = aes(x = long, y = lat, group = group), color = "gray") +
geom_point(data = pops_df, mapping = aes(x = lon, y = lat, col = ECS_SUBSEC, group = name)) +
scale_color_manual(values = ECS_cols) +
coord_fixed(ratio = 1.5) +
theme_bw()
ggplotly(MNsites)Circuitscape scripts were run on a cluster. See scripts in the scripts/samples resistance-ga circuitscape scripts for R
#Circuitscape rasters
precip_warm_q_rast <- stack("data/circuitscape_results/ind_varb_results/precip_warm_q_resist.tif")
proj4string(precip_warm_q_rast) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")
drainage_rast <- stack("data/circuitscape_results/ind_varb_results/drainage_resist.tif")
proj4string(drainage_rast) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")
imperv_rast <- stack("data/circuitscape_results/ind_varb_results/imperv_resist.tif")
proj4string(imperv_rast) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")
land_class_rast <- stack("data/circuitscape_results/ind_varb_results/land_class_resist.tif")
proj4string(land_class_rast) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")
mean_temp_warm_q_rast <- stack("data/circuitscape_results/ind_varb_results/mean_temp_warm_q_resist.tif")
proj4string(mean_temp_warm_q_rast) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")
min_temp_cold_m_rast <- stack("data/circuitscape_results/ind_varb_results/min_temp_cold_m_resist.tif")
proj4string(min_temp_cold_m_rast) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")
soil_comp_pc1_rast <- stack("data/circuitscape_results/ind_varb_results/soil_comp_pc1_resist.tif")
proj4string(soil_comp_pc1_rast) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")
soil_comp_pc2_rast <- stack("data/circuitscape_results/ind_varb_results/soil_comp_pc2_resist.tif")
proj4string(soil_comp_pc2_rast) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")
soil_min_pc1_rast <- stack("data/circuitscape_results/ind_varb_results/soil_min_pc1_resist.tif")
proj4string(soil_min_pc1_rast) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")
soil_min_pc2_rast <- stack("data/circuitscape_results/ind_varb_results/soil_min_pc2_resist.tif")
proj4string(soil_min_pc2_rast) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")precip_warm_q_p <- levelplot(precip_warm_q_rast, col.regions = rev(magma(256)), margin = FALSE, main = "Precip Warmest Q") +
spplot(states.lines.crop, col.regions = "slateblue1")
drainage_p <- levelplot(drainage_rast, col.regions = rev(magma(256)), margin = FALSE, main = "Drainage") +
spplot(states.lines.crop, col.regions = "slateblue1")
land_class_p <- levelplot(land_class_rast, col.regions = rev(magma(256)), margin = FALSE, main = "Land Classification") +
spplot(states.lines.crop, col.regions = "slateblue1")
mean_temp_warm_q_p <- levelplot(mean_temp_warm_q_rast, col.regions = rev(magma(256)), margin = FALSE, main = "Mean Temp of Warmest Q") +
spplot(states.lines.crop, col.regions = "slateblue1")
min_temp_cold_m_p <- levelplot(min_temp_cold_m_rast, col.regions = rev(magma(256)), margin = FALSE, main = "Min Temp of Coldest Month") +
spplot(states.lines.crop, col.regions = "slateblue1")
soil_comp_pc1_p <- levelplot(soil_comp_pc1_rast, col.regions = rev(magma(256)), margin = FALSE, main = "Soil Comp PC1") +
spplot(states.lines.crop, col.regions = "slateblue1")
soil_comp_pc2_p <- levelplot(soil_comp_pc2_rast, col.regions = rev(magma(256)), margin = FALSE, main = "Soil Comp PC2") +
spplot(states.lines.crop, col.regions = "slateblue1")
soil_min_pc1_p <- levelplot(soil_min_pc1_rast, col.regions = rev(magma(256)), margin = FALSE, main = "Soil Mineral PC1") +
spplot(states.lines.crop, col.regions = "slateblue1")
soil_min_pc2_p <- levelplot(soil_min_pc2_rast, col.regions = rev(magma(256)), margin = FALSE, main = "Soil Mineral PC2") +
spplot(states.lines.crop, col.regions = "slateblue1")
imperv_p <- levelplot(imperv_rast, col.regions = rev(magma(256)), margin = FALSE, main = "Impervious Surfaces") +
spplot(states.lines.crop, col.regions = "slateblue1")
paste("Precipitation in the warmest quarter: Resistance appears lower in the eastern half of Minnesota where it is wetter in the summer. The data transformation used was an inverse-reverse Ricker.")## [1] "Precipitation in the warmest quarter: Resistance appears lower in the eastern half of Minnesota where it is wetter in the summer. The data transformation used was an inverse-reverse Ricker."
precip_warm_q_ppaste("Mean temperature of the warmest quarter: Resistance appears lower in the northern half of Minnesota where it is cooler in the summer. The data transformation used was an inverse-reverse Ricker.")## [1] "Mean temperature of the warmest quarter: Resistance appears lower in the northern half of Minnesota where it is cooler in the summer. The data transformation used was an inverse-reverse Ricker."
mean_temp_warm_q_ppaste("Minimum temperature of the coldest month: Resistance appears lower in the southern half of Minnesota it is warmer in the winter. The data transformation used was an inverse-reverse Ricker.")## [1] "Minimum temperature of the coldest month: Resistance appears lower in the southern half of Minnesota it is warmer in the winter. The data transformation used was an inverse-reverse Ricker."
min_temp_cold_m_ppaste("Soil Composition PC1: Resistance is higher in the northern half of Minnesota where there are sandier soils and lower in the south with clay soils. The data transformation was a reverse Ricker.")## [1] "Soil Composition PC1: Resistance is higher in the northern half of Minnesota where there are sandier soils and lower in the south with clay soils. The data transformation was a reverse Ricker."
soil_comp_pc1_ppaste("Soil Composistion PC2: Resistance is realtively high for this entire layer with only a moderate drop in the central MN where there are areas of slightly higher PC2 values that indicate greater bulk volume density. The data transformation was a Ricker function.")## [1] "Soil Composistion PC2: Resistance is realtively high for this entire layer with only a moderate drop in the central MN where there are areas of slightly higher PC2 values that indicate greater bulk volume density. The data transformation was a Ricker function."
soil_comp_pc2_ppaste("Soil Mineral PC1: Resistance is generally low across the state and is slightly lower in eastern Minnesota where there is greater nitrogen. The data transformation was an inverse Ricker function.")## [1] "Soil Mineral PC1: Resistance is generally low across the state and is slightly lower in eastern Minnesota where there is greater nitrogen. The data transformation was an inverse Ricker function."
soil_min_pc1_ppaste("Soil Mineral PC2: Resistance is generally high across the state and is slightly lower in the Rainy river/Rainy Lake region where there is greater carbon content. The data transformation was an Ricker function.")## [1] "Soil Mineral PC2: Resistance is generally high across the state and is slightly lower in the Rainy river/Rainy Lake region where there is greater carbon content. The data transformation was an Ricker function."
soil_min_pc2_ppaste("Drainage Class: Resistance is especially high for the more poorly drained soils.")## [1] "Drainage Class: Resistance is especially high for the more poorly drained soils."
drainage_ppaste("Land Classification: Resistance is especially high for cropland/culitvated areas, hay/pastures, and grasslands.")## [1] "Land Classification: Resistance is especially high for cropland/culitvated areas, hay/pastures, and grasslands."
land_class_ppaste("Percentage impervious surfaces: Impervious surfaces had lower resistance.")## [1] "Percentage impervious surfaces: Impervious surfaces had lower resistance."
imperv_p#Genetic Chord Distance
tansy_gen_dist_mat <- as.dist(scale(as.matrix(read.csv("data/MN_tansy_gen_chord_dist_matrix.csv", header = T))))
#Spatial Distances
MN_geodist <- as.dist(fields::rdist.earth(xy_MN@coords))#Read in and scale resistance matrices then convert to 'dist' matrices for analysis
#pure resistance base on distance matrix
dist_mat <- as.dist(scale(as.matrix(read.csv("data/circuitscape_results/ind_varb_results/drainage/Distance_jlResistMat.csv", header = F))))
#Resistance matrices based on circuitscape optimizations
drainage_rast <- as.dist(scale(as.matrix(read.csv("data/circuitscape_results/ind_varb_results/drainage/drainage_bb_jlResistMat.csv", header = F))))
imperv_rast <- as.dist(scale(as.matrix(read.csv("data/circuitscape_results/ind_varb_results/imperv_surf/imperv_jlResistMat.csv", header = F))))
land_class_rast <- as.dist(scale(as.matrix(read.csv("data/circuitscape_results/ind_varb_results/Land_class_results/land_class_jlResistMat.csv", header = F))))
mean_temp_warm_q_rast <- as.dist(scale(as.matrix(read.csv("data/circuitscape_results/ind_varb_results/mean_temp_warm_q/mean_temp_warm_q_jlResistMat.csv", header = F))))
min_temp_cold_m_rast <- as.dist(scale(as.matrix(read.csv("data/circuitscape_results/ind_varb_results/min_temp_cold_m/min_temp_cold_m_jlResistMat.csv", header = F))))
soil_comp_pc1_rast <- as.dist(scale(as.matrix(read.csv("data/circuitscape_results/ind_varb_results/soil_comp_pc1/soil_comp_pca1_transpose_bb_jlResistMat.csv", header = F))))
soil_comp_pc2_rast <- as.dist(scale(as.matrix(read.csv("data/circuitscape_results/ind_varb_results/soil_comp_pc2/soil_comp_pca2_transpose_bb_jlResistMat.csv", header = F))))
soil_min_pc1_rast <- as.dist(scale(as.matrix(read.csv("data/circuitscape_results/ind_varb_results/soil_min_pc1/soil_min_pca1_transpose_bb_jlResistMat.csv", header = F))))
soil_min_pc2_rast <- as.dist(scale(as.matrix(read.csv("data/circuitscape_results/ind_varb_results/soil_min_pc2/soil_min_pca2_transpose_bb_jlResistMat.csv", header = F))))multisurf_run1_rast <- stack("data/circuitscape_results/multi_varb_results/multi-surface_results/clim + landclass + imperv + soil/MN_clipped_run1_multisurf_resist.tif")
proj4string(multisurf_run1_rast) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")
multisurf_run1_rast_scale <- scale(multisurf_run1_rast)
multisurf_run2_rast <- stack("data/circuitscape_results/multi_varb_results/multi-surface_results/clim + landclass + imperv + soil/MN_clipped_run2_multisurf_resist.tif")
proj4string(multisurf_run2_rast) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")
multisurf_run2_rast_scale <- scale(multisurf_run2_rast)
multisurf_run3_rast <- stack("data/circuitscape_results/multi_varb_results/multi-surface_results/clim + landclass + imperv + soil/MN_clipped_run3_multisurf_resist.tif")
proj4string(multisurf_run3_rast) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")
multisurf_run3_rast_scale <- scale(multisurf_run3_rast)
multisurf_avg_rast <- mean(multisurf_run1_rast_scale, multisurf_run2_rast_scale, multisurf_run3_rast_scale)multisurf_run1_p <- levelplot(multisurf_run1_rast_scale, col.regions = rev(magma(256)), margin = FALSE, main = "Multivariable Resistance Surface: Run 1") +
spplot(states.lines.crop, col.regions = "slateblue1")
multisurf_run2_p <- levelplot(multisurf_run2_rast_scale, col.regions = rev(magma(256)), margin = FALSE, main = "Multivariable Resistance Surfaces: Run 2") +
spplot(states.lines.crop, col.regions = "slateblue1")
multisurf_run3_p <- levelplot(multisurf_run3_rast_scale, col.regions = rev(magma(256)), margin = FALSE, main = "Multivariable Resistance Surfaces: Run 3") +
spplot(states.lines.crop, col.regions = "slateblue1")
multisurf_avg_p <- levelplot(multisurf_avg_rast, col.regions = rev(magma(256)), margin = FALSE, main = "Multivariable Resistance Surfaces: Average") +
spplot(states.lines.crop, col.regions = "slateblue1")
multisurf_run1_pmultisurf_run2_pmultisurf_run3_pmultisurf_avg_p_pts <- levelplot(multisurf_avg_rast, col.regions = rev(magma(256)), margin = FALSE, main = "Multivariable Resistance Surfaces: Average") + layer(sp.points(xy_MN, cex = 0.35, col = 1, pch = 19)) +
spplot(states.lines.crop, col.regions = "slateblue1")
multisurf_avg_p_pts